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1  INTRODUCTION 


The  social  and  economic  costs  associated  with  global  warming  will  be 
measured  in  terms  of  changes  in  the  frequency  and  intensity  of  extreme  events 
such  as  droughts,  floods,  hurricanes,  tidal  surges,  etc.  Small  changes  in  the 
mean  temperature  are  easily  adjusted  to,  but  successive  years  of  drought  or 
a  sequence  of  storms  are  much  less  easily  handled  even  by  advanced  societies. 
The  capability  to  forecast  the  frequency  and  intensity  of  extreme  events  in 
an  altered  atmosphere  poses  a  great  challenge. 

In  the  statistical  literature  the  large  or  small  values  assumed  by  a  ran¬ 
dom  variable  from  a  finite  set  of  measurements  are  termed  extreme  values. 
The  variable  of  interest  could  be  sea  level,  atmospheric  temperature,  pre¬ 
cipitation,  stream  flow,  etc.  The  largest  and  smallest  value  of  the  extremes 
are  often  of  most  interest.  However,  other  extremes  such  as  the  second  or 
third  largest  or  smallest  may  be  of  interest.  The  extreme  values  are  random 
variables.  For  example,  the  minimum  temperature  in  January  in  Bismarck, 
North  Dakota,  exhibits  random  variations  that  are  best  described  in  terms 
of  probabilities. 

The  application  of  extreme  value  statistics  is  uncommon  in  the  envi¬ 
ronmental  sciences  except  in  the  field  of  hydrology,  where  the  concepts  have 
been  used  in  designing  dams.  As  noted  in  Levine  et  al.  (1990),  extreme  value 
theory  has  not  received  much  attention  in  climate  studies  or  in  discussions  of 
greenhouse  warming  despite  the  obvious  importance  of  large  deviations  from 
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the  mean.  The  theory  and  applications  are  rarely  treated  in  books  on  prob¬ 
ability  and  statistics  and  almost  never  covered  in  university  courses.  If  there 
are  references  to  the  subject,  it  is  generally  to  Gumbel’s  (1958)  book  which 
has  long  been  out  of  print.  Cramer  and  Leadbetter  (1967)  and  Leadbetter, 
Lindgren  and  Rootzen  (1982)  provide  an  extensive  treatment  of  extremes  in 
stationary  sequences  but  focus  primarily  on  certain  of  the  difficult  mathe¬ 
matical  issues  related  to  extreme  values.  The  statistics  of  extremes  is  closely 
connected  to  order  statistics  which  provides  an  approach  to  distribution  free 
or  robust  statistics  (David,  1981).  The  present  paper  provides  a  review  of 
the  theory  with  a  number  of  applications  that  are  related  to  climate  change 
issues. 

The  applications  are  both  to  real  data  -  time  series  of  temperature 
derived  from  various  sources  -  and  to  results  obtained  from  much  simplified 
computer  models  of  climate.  The  latter  application  is  primarily  to  provide  an 
illustration  of  how  the  problem  of  extremes  can  be  approached,  rather  than 
to  give  definite  answers.  Some  of  the  examples  could  be  extended  to  data 
derived  from  large  scale  global  circulation  models  (GCMs)  of  the  atmosphere, 
but  since  we  do  not  have  access  to  GCMs  we  have  not  done  so. 

Sections  2,  3  and  4  serve  to  introduce  the  nomenclature  of  order  statis¬ 
tics  with  a  few  elementary  examples  of  the  calculation  of  the  probability  of 
extreme  climatic  events.  In  Section  5,  results  from  order  statistics  are  applied 
to  calculate  the  probability  that  the  observed  global  average  surface  air  tem¬ 
perature  records  contain  a  trend.  Unlike  a  previous  analysis  (Levine  et  al., 
1990),  no  assumption  is  made  as  to  the  underlying  statistical  distribution. 
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Since  extremes  are  rare  events  they  can  be  approximated  as  independent 
with  negligible  correlation.  In  accord  with  the  earlier  analysis  (Levine  et  al., 
1990),  order  statistics  indicate  very  high  odds  (105  —  106  to  1)  in  favor  of  the 
hypothesis  that  the  data  records  contain  a  linearly  increasing  trend  in  global 
average  surface  air  temperature. 

Section  6  presents  an  introductory  discussion  of  the  distribution  of  ex¬ 
tremes  with  no  restricting  assumptions  as  to  the  distribution  of  the  parent 
population.  While  general  results  can  be  obtained,  detailed  applications 
require  assumptions  as  to  the  underlying  distribution.  The  theory  is  spe¬ 
cialized  in  Section  7  to  normal  distributions.  An  application  to  the  Lorenz 
27-variable  model  of  climate  shows  that  the  temperature  values  calculated 
for  the  Lorenz  model  closely  approximate  normally  distributed  variates  both 
about  the  mean  and  in  the  tails  (extremes)  of  the  distribution  for  a  run  of 
1000  years. 

For  a  normally  distributed  random  variable,  the  range  (maximum  mi¬ 
nus  minimum)  increases  without  limit  as  the  number,  n,  of  variates  selected 
increases  (as  (In  n)1^2).  Section  7.2  provides  an  example  where  the  phys¬ 
ical  dimension  of  the  underlying  attractor  limits  the  growth  in  range  with 
time.  This  limit  is  not  observed  in  model  runs  or  in  observed  temperature 
records  (see  Sections  7.3  and  7.4).  Both  model  runs  and  observed  tempera¬ 
ture  records  from  which  a  trend  has  been  removed  closely  approximate  nor¬ 
mally  distributed  variates  about  the  mean  and  in  the  tails  of  the  distribution. 
These  results  raise  questions  as  to  the  predictability  of  climate. 

The  effect  of  a  shift  in  mean  on  the  frequency  of  extremes  is  discussed 
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in  Section  8  with  special  reference  to  possible  thresholds  for  damage  due  to 
climate  variability.  Section  9  provides  a  brief  summary  of  the  major  findings 
of  the  report. 
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2  RETURN  PERIOD 


Engineers  who  design  dams  are  interested  in  the  time  interval  between 
two  discharges  of  a  river,  each  of  which  is  greater  than  or  equal  to  a  given 
discharge.  The  probability  p  of  a  random  variate  X  having  a  value  greater 
than  x  is 

P  =  1  ~  F(x)  =  1  -  /  f(y)  dy 

J  —oo 

F{x)  =  Pr{  X  <  x ) 

where  F(x)  is  the  cumulative  probability  function  and  f(x)  is  the  probability 
density  function.  We  are  interested  in  the  number  of  independent  trials  k 
before  the  value  x  is  exceeded.  The  variable  &  is  an  integer  limited  on  the 
left  but  not  on  the  right,  since  the  value  of  x  may  never  be  surpassed.  The 
probability  that  x  is  first  exceeded  at  trial  k  is 

P(i-p)‘-‘  =  p«*-' 


since  the  event  failed  for  the  first  k-1  trials  (Bernoulli  trials).  The  mean  for 
k  is  simply 

The  return  period  is  defined  by 


T(x) 


1 

i  -F{xy 


If  an  event  has  probability  p,  on  the  average  we  have  to  make  1/p  trials  in 
order  for  the  event  to  happen  once. 
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The  standard  deviation  for  k  is 

a=V2(T2-  T)1/2 

which  for  a  long  return  period  becomes 

V?(T- i). 

The  smaller  the  probability  of  an  event,  the  larger  the  spread  of  the  distribu¬ 
tion.  This  general  property  lessens  the  usefulness  of  extremes  in  numerous 
situations. 

The  cumulative  probability  W(Jb)  that  the  event  happens  before  or  at 
the  kth  trial  is 

WW  =  i  -  (i  -  p)“- 

The  probability  that  the  event  will  happen  before  its  return  period  T  is 

W(T)  = 

a  l--  =  0.6231. 
e 

If  the  value  to  be  exceeded,  x,  is  large  and  the  return  period  long,  say 
greater  than  10,  then  the  cumulative  probability  W(£)  can  be  approximated 
by 

W(i)~l-exp(-^). 

In  hydrology  the  return  period  is  usually  measured  in  years.  For  exam¬ 
ple,  a  design  engineer  may  wish  to  have  a  high  probability  (0.95)  that  an 
event  does  not  happen  in  N  years.  The  return  period  T  for  the  design  is  then 
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approximately 


T  = 


N 

1  -  W 


so  that  in  order  to  have  a  probability  of  1  in  20  that  the  event  does  not 
happen  in  10  years  requires  the  design  to  have  a  much  larger  return  period 
of  200  years. 


These  simple  considerations  illustrate  certain  features  of  the  statistics  of 
extremes.  Under  the  assumption  that  the  events  are  independent,  results  can 
be  obtained  that  are  free  of  any  assumption  with  respect  to  the  underlying 
distribution.  The  assumption  of  independence  might  appear  to  be  highly 
restrictive.  For  dependent  variables,  their  correlation  as  a  function  of  time 
difference  or  lag  is  a  partial  measure  of  their  dependence.  Many  of  the  results 
derived  for  independent  variables  apply  to  dependent  variables  provided  that 
the  correlation  decreases  rapidly  enough  with  lag  (Leadbetter,  Lindgren  and 
Rootzen,  1983).  The  basic  reason  for  applicability  of  results  obtained  from 
the  theory  of  independent  events  to  correlated  series  is  that  extreme  events 
are  rare.  Rare  events  are  expected  to  be  separated  by  long  time  intervals 
so  that  the  effects  of  correlation  are  negligible.  If  the  correlation  function 
decays  at  least  as  fast  as  1  /In  n  where  n  is  the  time,  then  the  effects  of 
the  correlation  can  be  neglected  as  is  discussed  by  Leadbetter,  Lindgren  and 
Rootzen  (1983). 


7 


3  FREQUENCY  OF  EXCEEDANCES 


3.1  Bernoulli  Trials 


The  statistics  of  Bernoulli  trials  can  easily  be  applied  to  extreme  events. 
The  probability  that  during  the  next  n  years  the  maximum  July  temperature 
in  Washington,  D.C.,  will  exceed  105°F  exactly  k  times  is 

Pk  =  {  l  y-*( i-p)‘ 

where  p,  estimated  on  past  observations,  is  the  probability  that  the  maximum 
July  temperature  will  not  exceed  105°F  and  where 

V  k  '  (n  —  fc)!  k\ 

is  the  binomial  coefficient.  The  probability  of  exceeding  105°F  in  Washington 
is  about  a  1  in  20-  year  event,  p  =  0.95.  In  the  next  decade,  provided  that 
there  is  no  upward  or  downward  trend  in  temperatures,  the  probability  that 
the  maximum  July  temperature  will  not  exceed  105°F  is 

P0  =  (  ^  )(0.95)10  (0.05)°  =  0.6 

or  there  is  a  40  percent  chance  that  the  maximum  July  temperature  will 
exceed  105°F  in  the  period  1990-1999. 

The  probability  that  during  the  next  ten  years  there  will  be  fewer  than 
three  years  for  which  the  maximum  temperature  for  July  exceeds  105°F  can 
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also  be  calculated.  There  are  only  three  ways  that  the  number  of  exceedances 
of  105°F  cannot  be  greater  than  two  —  namely  when  the  number  is  zero,  one 
or  two.  The  probability  of  these  mutually  exclusive  events  is  the  sum  of  the 
probabilities  of  each 

Po  +  Px  +  Pi  =  (  q0  )  (0.95)10  (0.05)°  +  (  ^  )  (0.95)9  (0.05) 

+  (  10  )  (0.95)8  (0.05)2  =  0.9885. 

The  probability  of  the  maximum  July  temperature  in  Washington  in  the 
next  decade  exceeding  105°F  three  or  more  times  is  only  1  in  a  100.  If  in  the 
next  decade  the  temperature  does  exceed  105°F  in  July  three  or  more  times, 
the  hypothesis  that  there  is  no  upward  trend  in  the  temperature  should  be 
reexamined. 

The  number  of  exceedances  depends  on  the  interval  over  which  the  ob¬ 
servations  are  made.  From  the  results  of  the  above  example  we  can  calculate 
the  probability  that  over  the  next  century  every  decade  will  have  no  more 
than  two  Julys  in  which  the  maximum  temperature  exceeds  105°F.  With  a 
shift  in  time  scale  the  probability  is 

P  =  (  q0  )  (0.9885)10  (1  -  0.9885)°  =  0.89. 
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4  NON-PARAMETRIC  ORDER  STATIS¬ 
TICS 


The  need  to  estimate  the  probability  pm  of  the  mth  ranked  observation 
can  be  avoided  by  using  past  observations  to  determine  the  exceedances  in 
future  trials.  Let  n  denote  the  number  of  past  observations,  N  the  number 
of  future  observations,  m  the  rank  of  the  ordered  past  observation  and  k 
the  number  of  exceedances  of  the  mth  ranked  observation  in  N  future  trials. 
The  probability  of  observing  in  N  trials,  k  exceedances  of  the  mth  ranked 
observation  seen  in  n  trials  is: 


P(k;n, mtN)  =  (  *  )(  "  )mjf‘  pN-k  (1  -p)‘p”-m  (1  -p)m“Vp 


,  n  v  N  N 

< m  )( t )m 


since  the  integral  is  a  special  case  of  the  Beta  integral 

mm 

TV  +  *)  ‘ 


0(j,k)=  l'(\  ~yy~'  yk~1dy  = 
Jo 


For  the  largest  value,  m  =  1,  the  probability  of  k  exceedances  is 

nN\{N  +  n  -  k  -  1)} 


P(k;n,l,N)  = 


(N  -k)\  (N  +  n)  (N  +  n  —  1)!‘ 


The  probability  that  the  maximum  value  observed  in  n  trials  is  not  exceeded 
in  N  future  trials  is 


P(0;n,  l,N) 


n 

N  +  n 
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and  there  is  0.5  probability  that  the  maximum  is  not  exceeded  in  an  equal 
number,  N  =  n,  of  future  trials.  The  calculation  of  the  probabilities  for  the 
exceedance  of  the  largest  value  is  aided  by  the  recursion  relation 

P(k  +  1;  t.,  1,  JV)  =  ( ■?  |  *  ~  £  P(k-, n,  1  ,N). 


The  probability  that  all  values  in  N  future  trials  will  exceed  the  largest 
value  observed  in  the  first  n  trials  is  small.  For  the  case  N  =  n,  this  proba¬ 
bility  is 

which  for  N  =  6  is  0.00108. 


The  various  moments  of  the  distribution  function  can  be  obtained  from 
properties  of  the  hypergeometric  function  (Gumbel,  1954).  The  mean  number 
of  exceedances  is 

N 

E  (n,m,N)  =  k  P(k]n,m,N) 

fc=i 

=  mN/(n  +  1). 


Similarly,  the  variance  of  the  estimate  of  the  number  of  exceedances  is 

m(n  —  m  +  1)  N(N  +  n  +  1) 

(n  +  l)2(n  +  2) 

The  variance  increases  with  the  number,  N ,  of  future  triads  and  decreases 
with  the  number  of  past  observations  used  to  set  up  the  statistics.  The 
variance  for  the  number  of  exceedances  of  the  maximum  is 


<72  (n,m,  TV) 


<r2(n,l,JV) 


nN(N  +  n  +  1) 
(n  +  l)2  (n-f  2)' 
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The  variance  of  the  estimate  for  the  median,  m  =  n  odd,  is 

"  +  1  m  _  N(M  +  n  +  1) 

2  ’  ’  4(n  +  2) 

so  that  the  ratio  of  the  variance  of  the  estimated  number  of  exceedances  of 
the  maximum  to  the  number  of  exceedances  of  the  median  is 

<r2(n,  l,jV)  _  4n 

ai(N,^,N)  ~  (^TTT)2' 

In  one  sense,  the  extremes  are  more  reliable  than  the  median  since  the  vari¬ 
ance  of  the  number  of  exceedances  of  the  median  is  about  as  large  as  the 
variance  for  the  maximum.  In  the  limit  of  large  N  with  N  =  m  the  mean 
and  the  variance  take  on  simple  forms 

2?[fc]  ~  m 

and 

<r2(n,m,  N)  —  2m. 

The  mean  number  of  exceedances  over  the  mth  largest  value  is  equal  to 
the  rank  m  itself.  The  distribution  is  then  similar  to  a  Poisson  distribution 
for  integers.  However,  the  variance  is  twice  that  of  a  Poisson  distribution. 
The  probability  distribution  of  exceedances  in  this  case  is  independent  of  the 
distribution  of  the  random  variable  X.  It  does  depend  on  the  assumption  that 
the  observations  are  independent.  The  lack  of  dependence  on  the  distribution 
is  an  attractive  feature  of  the  theory,  since  the  distribution  may  be  known 
in  the  vicinity  of  the  median  but  there  will,  in  most  cases,  be  very  few 
observations  with  regard  to  extreme  values. 


<r2(n, 


5  APPLICATIONS  OF  ORDER  STATIS¬ 
TICS  TO  GLOBAL  AVERAGE  TEMPER¬ 
ATURE 


A  109-year  record  of  the  annual  global  average  temperature  of  the  at¬ 
mosphere  has  been  determined  by  Hansen  and  Lebedeff  (1987,  1988)  (see 
Figure  1).  The  data  have  been  normalized  to  zero  mean.  The  maximum 
temperature  for  the  100-year  period  1880-1979  is  0.2°C,  observed  in  1954. 
The  probability  of  one  or  more  exceedances  of  the  maximum  of  this  100- year 
record  in  the  nine- year  period  1980-1988  is 

P(>l;100,l,9)  =  l-^  =  0.083 

if  the  average  temperature  can  be  taken  as  an  independent  trial  and  there 
is  no  trend.  The  expected  number  of  exceedances  is  £(100, 1,9)  =  9/100  = 
0.089  with  a  variance  of  <r2(100, 1,9)  =  0.095.  In  fact  there  were  five  ex¬ 
ceedances:  1988,  1981,  1987,  1983  and  1980.  The  probability  under  the 
independence  assumption  of  five  exceedances  is 

P  (5;  100, 1,9)  =  1.03  x  10"6. 

The  observed  five  exceedances  strongly  contradict  the  assumption  of  inde¬ 
pendence  and  provide  support  for  the  hypothesis  that  there  is  an  increasing 
trend  in  the  annual  global  average  surface  temperature. 

If  a  linear  trend  with  a  least  square  slopes  of  0.55°  per  century  is  removed 
the  maximum  residual  in  the  first  100  years  is  0.28°C  in  1926.  There  are  no 
exceedances  for  the  years  1980-1988  in  the  residual  as  would  be  expected  from 
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Table  1 


Expected  Number  of  Exceedances  in  Six-Year  Runs  of  Global 
Annual  Average  Temperature  and  Observed  Exceedances  in  17 

Six- Year  Runs 


Number  of 
Exceedances  k 
of  Maximum 
Value 

Probability 
of  k 

Expected 
Number  of 
Exceedances 

Observed  Number 
of  Exceedances 

Observed  Number 
of  Exceedances 
in  Residual 
after  Removal 
of  Linear 

Trend 

0 

0.5 

8.5 

0 

2 

1 

0.273 

4.64 

1 

3 

2 

0.136 

2.31 

1 

3 

3 

0.061 

1.03 

0 

4 

4 

0.027 

0.46 

2 

2 

5 

0.006 

0.11 

1 

1 

6 

0.001 

0.02 

12 

2 

the  probability  distribution  for  the  maximum.  In  terms  of  odds,  the  ratio  of 
the  probability  of  the  observations  containing  a  trend  to  the  probability  that 
the  temperatures  are  independent,  identically  distributed  random  variables 
is 


0.917 

1.03  x  10-6 


=  8.9  x  105. 


As  a  further  example  of  the  use  of  the  probability  distribution,  the  108- 
year  record,  1980-1987,  is  divided  into  18  consecutive  six-year  segments.  The 
maximum  value  in  the  first  six-year  period,  1880-1885  is  -  0.37°C.  The  num¬ 
ber  of  exceedances  in  each  of  the  subsequent  six-year  segments  is  tabulated 
in  Table  1  along  with  the  number  obtained  for  the  probability  distribution. 
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The  expected  number  of  exceedances  in  any  six-year  period  is 

E  (6, 1,6)  =  6/7  =  0.86 
with  a  standard  deviation  of 

a  =  1.09. 

Clearly  the  number  of  six-year  runs  with  all  six  values  in  excess  of  the  maxi¬ 
mum  for  the  first  six-year  segment  is  much  larger,  16,  than  would  be  expected 
from  a  sequence  of  independent  trials.  The  probability  of  obtaining  12  seg¬ 
ments  out  of  17  in  which  all  the  values  are  in  excess  is 

171 

— ^  (0.00108)12  (1  -  0.00108)5  =  1.5  x  1(T32. 

12!  5! 

Table  1  also  tabulates  the  observed  distribution  of  exceedance  for  the  resid¬ 
uals  after  a  linear  trend  has  been  removed.  The  observed  distribution  still 
differs  from  the  calculated  one.  However,  the  probability  of  obtaining  two 
segments  in  which  all  six  values  axe  greater  than  the  maximum  observed  in 
the  first  six-years  while  small,  1.5  x  104,  is  much  larger  than  the  probability 
of  obtaining  12  such  sequences. 

Levine  et.  al  (1990)  calculated  that  the  odds  for  a  linear  trend  in  the 
data  were  1.7  x  105  to  one.  In  that  analysis,  the  correlation  between  annual 
temperature  was  explicitly  taken  into  account  through  the  use  of  the  full 
covariance  matrix  in  the  least  squares  analysis  but  the  analysis  was  based 
on  the  assumption  that  the  global  average  temperatures  were  normally  dis¬ 
tributed.  The  analysis  presented  above  supports  the  hypothesis  of  a  trend  by 
similar  odds,  but  is  based  on  the  assumption  of  independent  trials  without 
any  assumption  on  the  underlying  statistical  distribution. 
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6  DISTRIBUTION  OF  EXTREMES 


The  above  analysis  provides  information  with  respect  to  the  number  of 
observations  that  would  exceed  the  mth  largest  value  of  n  initial  observations. 
The  analysis  provides  no  information  on  the  value  that  would  exceed  the 
observed  extremes.  Such  information  requires  that  the  underlying  probability 
distribution,  F(x),  be  known.  As  illustrated  above,  some  of  the  analysis  of 
the  frequency  of  exceedances  can  be  carried  out  without  assumptions  with 
respect  to  the  underlying  distribution  of  the  random  variable.  The  exact 
form  of  the  distribution  is  generally  not  known,  but  an  asymptotic  theory 
that  covers  a  wide  variety  of  initial  distribution  has  been  developed  (Cramer, 
1946,  Gumbel,  1958). 


As  before,  let  k  denote  the  kth  value  from  the  top  of  a  sample  of  n 
observations.  The  probability  density  g(x]  n,  k)  for  the  fcth  value  is  equal  to 
the  probability  that  among  n  samples  n  —  k  are  less  than  x,  and  k  —  1  are 
greater  than  x  while  the  remaining  value  falls  between  x  and  x  +  dx 

«(*;»■*)  =  »(  “  "  }  )(F(x r*  (1  -  F(x))‘-']/(x)<fx 

where  F(x)  is  the  probability  distribution  function  and  /(x)  is  the  density 
distribution  function.  For  the  maximum  value  k  =  1  and 


0(x;n,l)  =  (F(x)]"  dx. 


For  a  population  having  a  normal  distribution  with  zero  mean  and  unit  vari¬ 
ance 


F(x) 


1 


e 


~t7'2  dt. 
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The  dependence  of  the  maximum  value  on  the  number  of  observations  from 
a  series  of  independent  trials  drawn  from  a  normal  population  is  illustrated 
in  Figure  2.  For  a  sample  of  ten,  the  maximum  has  a  0.05  probability  of 
exceeding  2.57  standard  deviations;  for  a  sample  of  1000,  the  maximum  will 
exceed  3.89  standard  deviations  one  in  twenty  times. 
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Figure  2.  Variations  of  the  probability  of  the  maximum  of  n  independent  trials  with  n. 
The  distribution  has  zero  mean  and  unit  variance. 
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7  ASYMPTOTIC  DISTRIBUTION  OF  EX¬ 
TREMES  FOR  A  NORMAL  PROCESS 


The  general  asymptotic  theory  of  the  distribution  of  extremes  deals  with 
the  conditions  under  which 


P[an(Mn  —  «w)  <  *] 

approaches  a  limiting  distribution  G(x)  where  an  and  bn  are  suitable  normal¬ 
izing  constants  and 

Mn  max[^fj ,  .Afj,  •••  ,  Xn\  • 

The  theory  states  that  every  distribution  G,  has,  up  to  scale  and  location 
changes,  one  of  the  following  parametric  forms  commonly  called  the  Three 
Extreme  Value  distributions: 

Type  I:  oo  <  z  <  oo  G{x)  =  exp(— e~x) 

Type  II:  z  >  0  G(x)  =  exp(— z~°);  a  >  0 

z  <  0  G(x)  =  0 

Type  III:  z  <  0  G(x)  =  exp(— (— z)°);  o  >  0 

z  >  0  =1 

(Leadbetter,  Lindgren  and  Rootzen,  1983).  The  normal  distribution,  unlim¬ 
ited  on  the  left  and  right,  falls  into  the  Type  I  domain  of  attraction.  The 
density  distribution  for  the  fcth  ordered  value  becomes 

exp(—  xk  —  e-r) 

which  is  illustrated  in  Figure  3  for  k  =  1 , 2, 3. 

For  a  normal  distribution  with  zero  mean  and  unit  variance,  if  z  is  the 
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-2-10123456 

X 

Figure  3.  The  probability  density  function  for  the  double  exponential. 
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Jbth  value  from  the  top,  then  £  is 


n 

7%  J*  e 


-t2!2dt. 


The  normalizing  constants  can  be  obtained  by  a  complex  calculation. 
The  solution  for  x  when  n  is  large  can  be  obtained  by  partial  integration 
leading  to 

=  -  e +0(\)e~l3/2. 
n  x  x 2 

With  bounded  £,  x  is  given  by 

In  In  n  +  In  Air  In  £  1  . 

x  =  (2  n  n)  -  2(2  In  n )*/*  "  (2  In  n)1/2  +  °(7^‘ 

For  the  general  form  of  the  normal  distribution  with  arbitrary  mean  m  and 
variance  <r2,  the  form  for  M„  =  x  is  asymptotically 

_  ,  ,1/7  In  In  n  + In  Air  av 

x  =  m  +  <7(2  In  n)I/2  -  <7-^^ - rw^-  + 


2(2  Inn)1'2  (2/nn)'/J 

where  v  is  a  random  variable  having  the  frequency  function 

exp  (—xk  —  e~r). 

The  expression  for  the  kth  value  from  the  bottom  is  given  by  a  similar  ex¬ 
pression  with  opposite  sign  for  all  terms  but  the  mean,  m. 


The  expected  value  for  the  kth  value  from  the  top  is 

+  In  Air  +  2(Si  —  C) 


E\x\  n,  k)  =  m  +  a  ^2  In  n)1^2  —  - 


(2  In  n)1/2 


+  0( 


with  variance  D 2 


D2{x'n'k)  =  (y-Si)+0(/^;) 
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where  C  is  Euler’s  constant,  £7=0.57722..., 


„  1  1 
S,-I+2+ 

^  +  ^  + 


+ 


1 


k~l 

1 


Similar  expressions  can  be  calculated  for  the  difference  between  the  kth  max¬ 
imum  and  fcth  minimum  values  of  the  n  independent  trials  drawn  from  a 
normal  population. 


E[x  -  y;  n,k]  =  cr  —  — 


In  Inn  —  In  4x  —  2(5i  —  C) 
(2  In  n)1/2 


+  0( 


inn*) 


0‘[x-y,n,k)  =  ^—(’?r  -  5j)  +  Ot-j-l— )• 
Inn  6  In 2  n 


The  dependence  of  the  range  on  the  number  of  observations  is  indicated 
in  Figure  4.  The  three  standard  deviation  bounds  are  also  shown.  For  a 
normal  population  of  a  1000  variables  there  is  about  1  in  a  100  probability 
that  the  maximum  range  will  exceed  eight  standard  deviations.  The  expected 
value  for  the  maximum  of  normal  population  is  shown  in  Figure  5. 


7.1  Applications  of  Asymptotic  Theory  of  Extremes 
to  a  Climate  Model 


The  variation  of  global  average  annual  temperature  depicted  in  Figure  1 
shows  irregular  behavior  with  a  non-vanishing  correlation  function  for  lags 
up  to  the  length  of  the  record  and  continuous  power  spectrum  with  abun¬ 
dant  energy  at  low  frequencies  (Levine  et.al,  1990).  The  irregular  behavior, 
characteristic  of  chaotic  systems,  also  shows  clearly  in  complex  models  of  the 
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Expected  Value  of  Maximum  Range 
(Units  =  Standard  Deviations) 


Figure  4.  Variations  of  the  maximum  expected  range  for  a  random  variate  drawn  from  a 
normal  population  as  a  function  of  the  number  of  observations.  The 
distribution  has  zero  mean  and  unit  variance.  The  three  standard  deviations 
from  the  expected  value  are  also  shown. 
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Expected  Value  of  Maximum  Value 
(Unit  =  Standard  Deviation) 


0  200  400  600  800  1000 


Number  of  Observations 


Figure  5.  Variation  in  expected  maximum  value  of  an  observation  drawn  from  a  normal 
population  with  zero  mean  and  unit  variance  as  a  function  of  the  number  of 
observations.  The  three  standard  deviations  from  the  expected  value  are  also 
shown. 
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atmosphere  such  as  General  Circulation  Models,  which  contain  about  100,000 
coupled  nonlinear  equations,  and  also  in  the  simplest  models.  Lorenz  (1984a) 
has  developed  a  low-order  model  in  27  variables  that  has  been  analyzed  in 
some  detail  by  Levine  et  al.,(1990).  The  physical  variables  include  the  mix¬ 
ing  ratios  of  water  vapor  and  liquid  water  as  well  as  temperature,  pressure 
and  the  usual  dynamic  variables.  In  the  Lorenz  model,  the  atmosphere  and 
ocean  exchange  heat  and  water  through  evaporation  and  precipitation.  The 
model  also  produces  clouds,  which  reflect  incoming  solar  radiation,  while 
both  phases  of  water  absorb  and  re-emit  infrared  radiation. 

For  a  1000- year  run  for  the  global  average  annual  temperature,  the  ob¬ 
served  standard  deviation  was  0.3583° C.  The  data  scaled  to  unit  variance  and 
zero  mean  are  shown  in  Figure  6.  The  observed  maximum  value  is  3.407.  The 
probability  of  obtaining  this  value  with  1000  values  is  0.712,  or  a  probability 
of  0.29  that  3.4  would  be  exceeded  in  a  1000-year  run  of  independent  trials 
drawn  from  a  normal  population.  The  distribution  of  temperature  deviations 
about  the  mean  approximate  a  normal  distribution  as  illustrated  in  Figure  7. 
The  variation  of  range  with  the  number  of  observations  is  shown  in  Figure  8. 
The  variation  closely  follows  the  expected  value  of  the  range  drawn  from  a 
normal  population. 

7.2  Physical  Dimensions  of  an  Attractor 


For  a  truly  random,  non-deterministic  process,  such  as  a  variate  drawn 
from  a  distribution  unlimited  on  the  left  and  right,  the  range  will  grow  with 
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Temperature  ( °C) 
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Figure  6.  Variation  of  global-mean,  annual-mean,  sea-level  air  temperature  for  1000  years 
in  a  numerical  solution  to  the  Lorenz  27-variable  model  of  the  atmosphere. 

Data  are  rated  to  unit  variance  and  zero  mean. 
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Figure  7.  Distribution  of  global-mean,  annual  average  temperatures  for  1000  years  of  the 
Lorenz  model. 
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Time  (Years) 


Figure  8.  The  variation  with  time  of  the  range  of  global-mean,  annual  average 

temperature  for  the  Lorenz  model.  The  dotted  curve  gives  the  expected  value 
for  the  range  for  samples  drawn  from  a  normal  population. 
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the  number  of  variates  selected  at  a  rate  approximately 

2(2  In  n)1^2. 

For  a  time  series  n,  the  number  of  variates  equals  the  number  of  time  units. 
In  a  deterministic  system  governed  by  a  set  of  coupled  nonlinear  differen¬ 
tial  equations,  one  would  expect  that  the  embodied  conservation  equations 
would  set  limits  to  the  maximum  growth  of  any  state  variable.  In  terms  of  the 
language  of  nonlinear  dynamical  systems,  the  governing  attractor  has  finite 
physical  dimensions.  In  this  context  we  are  speaking  of  the  physical  dimen¬ 
sions  of  the  attractor  as  contrasted  to  the  usual  embedding  dimensions  which 
constitute  a  measure  of  the  denseness  of  the  set  making  up  the  attractor. 

In  the  limit  of  large  times,  the  range  of  random  variables  grows  without 
bound  but  very  slowly  as  (In  n)1^2.  For  deterministic  systems,  the  growth,  as 
noted  above,  should  be  limited  at  a  large  enough  n.  For  the  Lorenz  system, 
it  is  apparent  that  the  limit  has  not  been  reached  in  1000  years  (see  Figure  8). 
A  general  problem  is  that  of  determining  the  limits  to  the  range  of  a 
particular  variable  and  the  time  required  to  explore  the  far  reaches  of  the 
attractor.  It  would  appear  that  these  are  difficult  issues,  particularly  for  a 
system  as  complex  as  the  27-variable  Lorenz  system. 

In  order  to  further  examine  the  time  required  for  the  Lorenz  27-variable 
model  to  explore  the  outer  regions  of  the  attractor,  the  model  was  run  for 
10,000  years.  The  results  are  shown  in  Figure  9.  It  would  appear  that 
the  times  to  attain  the  outer  reaches  of  the  attractor  are  on  the  order  of 
5,000  years.  To  distinguish  between  a  random  process  with  the  events  drawn 
from  a  normal  population  and  a  deterministic  Lorenz  system  by  examining 


33 


Maximum  -  Minimum 


Globa)  Averaged  Tempera, u 


Ran°e,0f  Gaussian  Distribution 


i - 


Sea  Surface  Temperature 


<000  ~  '  rr — 

6000 

■firne  in  Years 


%«*  9-  The  variation  with  , 

'm'n  frora  ■  — -  Population, *Xpect€d  —  of 


100000 


34 


the  tails  of  the  distribution  one  would  need  to  select  on  the  order  of  104 
values.  This  observation  suggests  that  before  asymptotic  behavior  is  reached 
by  a  dynamical  system,  transients  of  considerable  duration  can  occur.  Also 
shown  in  Figure  9  is  the  variation  of  the  range  of  sea  surface  temperature 
and  dew  point  temperature  as  calculated  in  the  model.  The  tails  of  trails  of 
the  distribution  of  the  parameters  differ  from  that  of  a  normal  population, 
but  very  long  times,  on  the  order  of  104  years,  are  required  to  reach  the 
outer  edges  of  the  attractor.  There  is,  of  course,  no  guarantee  that  in  any 
single  integration  the  system  will  reach  the  boundary  of  the  attractor.  A 
further  consideration  is  the  random  roundoff  noise  generated  in  the  computer 
integration.  This  noise  will  cause  the  calculated  range  to  mimic  a  random 
variable. 

In  a  simple  system,  again  one  proposed  by  Lorenz  (1984b),  it  does  ap¬ 
pear  possible  to  obtain  bounds  on  the  maximum  dimensions  of  the  attractor. 
Lorenz  (1984b)  labels  the  following  set  of  three  coupled  nonlinear  equations 
as  the  simplest  possible  general  circulation  model: 

=  -Y2-Z2-aX  +  aF 
=  XY  -  bXZ  -  Y  +  G 
=  bXY  +  XZ  -  Z. 

The  variable  X  represents  the  strength  of  the  large-scale  westerly-wind  cur¬ 
rent  or  the  geostrophically  equivalent  equatorial-pole  temperature  gradient, 
while  Y  and  Z  are  the  strengths  of  the  cosine  and  sine  phases  of  a  chain  of 
superposed  waves  (Lorenz,  1984b).  The  quadratic  terms  containing  b  rep¬ 
resent  the  translation  of  the  waves  by  the  western  current.  The  remaining 


dX 

dt 

dY 

dt 

dZ 

dt 
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quadratic  terms  represent  a  continual  transfer  of  energy,  except  when  X  is 
negative,  from  the  westerly  flow  to  the  waves. 


The  principal  external  driving  force,  the  contrast  between  equatorial 
and  pole  heating,  acts  directly  on  the  zonal  flow  X  and  is  represented  by  a 
F.  A  secondary  driving  force,  dependent  on  the  difference  between  oceans 
and  continents,  enters  the  equations  through  G. 

A  section  of  the  attractor  for  the  Lorenz  84  model  is  displayed  in  Figure  10, 
which  shows  the  range  for  the  Z  variable  projected  onto  a  plane  including 
the  Z  axis  and  at  an  angle  0  to  the  X  axis.  The  motion  is  centered  about 
Z  =  0  but  otherwise  offset  in  X  and  Y . 

We  treat  X ,  Y  and  Z,  as  coordinates  in  three-dimensional  phase  space. 
A  rate  of  change  of  energy  equation  follows  from  the  governing  set  by  multi¬ 
plying  each  by  2X,2Y,  and  2 Z,  respectively  and  summing 

jdJ 

—  =  -[«( 2X  -  F)2  +  (2 Y  -  G)2  +  (2Z)2  -  ( aF 2  +  G2)]/2 

where  R 2  =  X2  +  Y2  +  Z2  can  be  interpreted  as  the  total  energy  -  the  sum 
of  the  kinetic,  potential  and  internal  energies.  The  right-hand  side  of  the 
equation  for  ~  vanishes  on  an  ellipsoid 

a(X  -  F/2)2  +  {Y-  G/2)2  +  Z2  =  . 

4 

The  rate  of  change  of  energy  outside  the  ellipsoid  is  negative.  This  implies 
that  for  any  sphere  centered  on  the  origin  in  state  space  enclosing  the  el¬ 
lipsoid,  orbits  passing  through  points  outside  the  sphere  will  eventually  pass 
through  the  sphere  and  remain  inside.  The  ellipsoid  is  centered  at  (F/2,  G/2, 
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21 


Figure  10.  A  section  of  the  attractor  for  the  Lorenz  3-variable  model  at  a  plane  containing 
the  z  axis  and  at  an  angle  of  117°  to  the  x  axis.  The  parameters  for  the  model 
are  given  in  Table  2. 
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Table  2 


Maxima  and  Range  for  Variables  Defined  in  Lorenz 
Lowest  Order  Central  Circulation  Model 


(F=8.0, 

G=1.0,  a= 

0.25,  b=4) 

Variable 

Maximum 

Minimum 

X 

12.24 

-0.12 

Y 

3.06 

-1.06 

Z 

2.06 

-2.06 

0)  with  the  corresponding  major  axes  (s£2±2i)i/ (°.£?+G2  j1/2). 

The  maxima  and  minima  for  the  state  variables  are  given  in  Table  2  for 
initial  conditions  lying  within  the  ellipsoid  and  for  specific  values  of  the  con¬ 
stants  F,  G  and  a.  For  initial  conditions  lying  within  the  ellipsoid  defined  by 
^  =  0,  all  orbits  would  be  expected  to  remain  within  the  ellipsoid,  except 
for  small  inertial  excursions,  in  the  absence  of  noise.  The  ellipsoid  provides 
the  limits  for  values  that  can  be  taken  on  by  the  state  variables.  The  con¬ 
stant  b  does  not  enter  into  the  equation  for  the  critical  energy  surface  since 
b  defines  the  strength  by  which  the  waves  are  carried  by  the  zonal  westerly 
current  and  does  not  alter  the  energy  of  the  system. 

The  above  analysis  provides  bounds  on  the  maxima,  minima  and  range 
attained  by  the  state  variables  but  provides  no  estimate  of  the  average  time 
required  for  a  state  variable  to  reach  the  boundary  of  the  attractor.  The  only 
time  scale  within  the  systems  of  equations  is  determined  by  the  constant  a, 
the  inverse  of  the  decay  time  for  the  zonal  current.  Lorenz  (1984b)  assumes 
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a  decay  time  for  the  zonal  current  of  20  days  by  selecting  a=0.25  and  taking 
5  days  as  the  unit  time. 

A  numerical  trial  employing  the  constants  shown  in  Table  2,  taking  the 
time  unit  to  be  6  hours,  averaging  over  2.5  days  and  integrating  from  34.25 
years  (5,000  points)  yielded  a  variation  in  the  X  coordinate  shown  in  Figure  1 1 . 
The  power  spectrum  and  autocorrelation  function  for  X ,  equivalent  to 
the  equator-pole  temperature  difference,  are  shown  in  Figures  12  and  13. 

The  power  spectrum  displays  much  more  structure  than  the  Lorenz  27- 
variable  model  even  though  there  is  only  one  time  scale  that  enters  through 
the  constant  a.  The  distribution  of  values  of  X  departs  markedly  from  that 
a  normal  distribution  (see  Figure  14)  in  contrast  to  the  Lorenz  27- variable 
model. 

The  maximum  range  observed  in  the  data  for  X  is  2.98,  well  within  the 
limits  shown  in  Table  2.  In  fact,  in  this  particular  run,  the  calculated  orbits 
did  not  explore  the  outer  edges  of  the  attractor  for  the  state  variable  X.  The 
variation  of  the  maximum  range  with  time  for  X  normalized  to  zero  mean 
and  unit  variance  is  shown  in  Figure  15.  Unlike  the  the  Lorenz  27-variable 
model,  the  maximum  range  does  not  continue  to  climb  with  time  but  levels 
of  well  within  the  limits  set  by  the  critical  energy  surface. 

The  Y  and  Z  state  variables  show  similar  behavior.  The  variation  of  Z 
with  time  is  shown  in  Figure  16.  The  histogram  (Figure  17)  for  the  Z  vari¬ 
able  shows  large  deviation  from  a  normal  population.  The  range  levels  off 
with  the  maximum  observed  range  of  4.17  (see  Figure  18)  slightly  exceeding 
the  maximum  range  expected  for  the  Z  variable  of  4.12  (see  Table  2).  An 
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Time  (unit  =  2.5  days) 


Figure  11. 


Variation  of  the  Lorenz  3-variable  model  X  coordinate  (X  is  the  equatorial-pole 
temperature  gradient)  with  time  for  12,500  days. 


40 


Figure  13.  Autocorrelation  junction  for  state  variable  X  in  Lorenz  3-variable  model.  The  data 
are  normalized  to  zero  mean  and  unit  variance. 
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Figure  14.  The  distribution  of  value  for  the  X  state  variable  in  the  Lorenz  3-variable  model. 
The  data  have  been  normalized  to  zero  mean  and  unit  variance. 
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Figure  15.  Variation  of  the  range  of  the  state  variable  X  in  the  Lorenz  3-variable  model. 

The  data  are  normalized  to  zero  mean  and  unit  variance.  The  dotted  line  is 
the  range  expected  for  values  chosen  from  a  normal  population  with  unit 
variance  and  zero  mean. 
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Figure  16.  Variation  of  the  Z  state  variable  with  time  for  the  Lorenz  3-variable  model. 
The  unit  of  time  is  2.5  days. 


Figure  17.  Histogram  for  the  Z  state  variable  in  the  Lorenz  1984  model. 
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Figure  18.  Range  for  the  Z  state  variable  in  the  Lorenz  1984  model.  Data  have  been 
normalized  to  unit  variance  and  zem  mean.  The  curve  represents  expected 
range  for  a  normally  distributed  var.able. 
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observation  of  a  value  for  the  range  slightly  in  excess  of  the  maximum  pre¬ 
dicted  range  is  not  unexpected,  since  inertia  will  carry  the  the  orbit  outward 
before  the  inward  pull  of  the  attractor  brings  the  orbit  back  within  the  limits 
of  the  critical  energy  surface. 

The  above  example  illustrates  that  it  is  possible  to  predict  the  maximum 
range  of  variables  in  a  complex  nonlinear  system  provided  the  equations  are 
known.  The  limitation  to  the  range  of  values  which  variables  can  assume 
likely  comes  from  the  conservation  equations  governing  the  motion.  The 
time  scale  over  which  the  maximum  value  may  be  attained  is  not  determined 
by  the  above  analysis. 

7.3  Extremes  for  a  Climate  Model  with  a  Linear  In¬ 
crease  in  Temperature 


Levine  etal.,  (1990)  describe  a  version  of  the  Lorenz  27-variable  model 
in  which  the  solar  insolation  is  increased  linearly  with  time  leading  to  an 
increase  of  global  temperature  of  10°  over  a  1000-year  period.  The  resulting 
temperature  record  is  shown  in  Figure  19.  The  calculated  variance  is  9.24, 
which  can  be  broken  down  into  the  part  due  to  the  linear  increase  and  a2, 
the  variance  about  the  linearly  increasing  mean. 


N2b2 

12 

0.235 


+  o 


2 


where  N  is  the  number  of  years  and  b  is  the  rate  of  temperature  increase. 
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Figure  19.  Global  average  temperature  for  Lorenz  2  7- variable  model  of  the  atmosphere 
in  which  the  parameter  corresponding  to  the  insolation  is  increased  linearly 
in  time  to  provide  an  approximately  10°C  increase  in  temperature  over 
1000  years. 
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The  residual  temperature  curve  after  removing  a  trend  by  least  squares 
is  shown  in  Figure  20  with  a  distribution  of  the  residuals  approximating  a 
nornr'!  distribution  (see  Figure  21).  The  maximum  range  after  1000  years 
for  the  residuals  approximates  the  range  expected  by  drawing  variates  from  a 
normal  population  (see  Figure  22).  The  removal  of  a  trend  before  analyzing 
for  extremes  is  a  useful  way  of  examining  records  of  global  average,  surface  air 
temperature  records.  There  is  a  high  probability  that  these  records  contain 
a  linear  trend  (Levine  et  al.  1990). 

7.4  Distribution  of  Extremes  for  Global  Average  Tem¬ 
perature  Record 


Our  final  example  of  an  application  of  extreme  value  theory  is  to  an  an¬ 
nual  average,  global  surface  air  temperature  record.  Three  research  groups 
(Jones,  1988;  Hansen  and  Lebedeff,  1987,  1988;  and  Vinnikow  et  al.,  1990) 
have  produced  similar  analyses  of  hemispheric  surface  temperature  variations 
from  somewhat  differing  initial  data  sets.  The  longest  of  the  data  sets  is  that 
of  Jones  (1988)  and  we  will  use  it  in  the  analysis.  An  analysis  of  the  Hansen- 
Lebedeff  (1987,  1988)  series  (see  Figure  1)  produced  results  similar  to  those 
obtained  from  the  Jones  record.  The  Jones  global  average  temperature  vari¬ 
ation  is  shown  in  Figure  23,  which  can  be  compared  with  Figure  1.  Both 
figures  show  a  relatively  cool  late  19th  century  followed  by  warming  inter¬ 
rupted  by  cooling  in  the  1940-70  period  (see  Levine  et  al.,  1990).  The  power 
spectrum  for  the  Jones  (1988)  surface  air  temperature  records  is  shown  in 
Figure  24.  The  distribution  of  values  for  the  Jones  record  is  shown  in  Figure  25. 
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Figure  21.  Distribution  of  values  of  the  temperature  for  record  shown  in  Figure  18. 
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Range  in  Temperature  (CC) 


Figure  22.  Range  for  the  residuals  for  the  Lorenz  27-variable  system  shown  in  Figure  20. 

The  data  have  been  normalized  to  zero  mean  and  unit  variance.  The  curve  is 
the  expected  range  for  a  normal  population. 
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Figure  23.  The  variation  of  annual  average  global  surface  air  temperature  according  to 
Jones  (1988). 
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Figure  24.  Power  spectrum  for  the  global  annual  average  surface  air  temperature  record 
shown  in  Figure  23.  The  data  have  been  normalized  to  zero  mean  and  unit 
variance. 
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The  distribution  differs  from  a  normal  distribution  by  being  flatter.  This 
is  consistent  with  the  observed  trend  of  increasing  temperature  displayed  in 
Figure  23  (a  straight  line  would  exhibit  a  uniform  distribution  of  values). 
The  range  for  the  Jones  record  (see  Figure  23)  is  shown  in  Figure  26.  The 
deviations  again  are  in  the  direction  expected  by  a  record  containing  a  linear 
trend. 

If  a  trend  determined  by  least  squares  (slope  =  0.28°C/century)  is  re¬ 
moved  from  the  record,  the  resulting  temperature  residuals  exhibit  a  his¬ 
togram  that  closely  approximates  a  normal  distribution  (see  Figure  27).  The 
tails  of  the  distribution  also  approximate  those  of  a  normal  variate,  as  is  illus¬ 
trated  in  Figure  28,  which  shows  the  variation  of  the  range  with  time  for  the 
record  exhibited  in  Figure  23  from  which  a  linear  trend  has  been  removed. 

The  above  analysis  leads  to  the  conclusion  that  once  the  linear  trend 
is  removed  from  the  observed  record,  the  distribution  of  the  residuals  both 
about  the  mean  and  in  the  tails  approximates  a  normal  distribution.  In  the 
Lorenz  27-variable  system  we  observed  chaotic  behavior  which  led  to  a  ran¬ 
dom  (normally  distributed)  variation  even  over  time  scales  of  thousands  of 
years.  The  observed  data,  even  when  temporally  and  spatially  averaged,  ex¬ 
hibit  similar  behavior.  These  observations  suggest  the  possibility  of  chaotic 
noise  at  low  frequency  with  important  implications  for  predictability.  In  par¬ 
ticular,  records  of  less  than  1000  years  may  show  statistics  indistinguishable 
from  a  random  series  (once  trends  assumed  to  be  deterministic  have  been 
removed).  Prediction,  outside  of  the  trend,  may  thus  be  impossible  except  in 
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Figure  25.  Histogram  for  global  annual  average  surface  air  temperature  record  shown  in 
Figure  23. 
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Figure  26.  Variation  of  range  for  the  global  annual  average  surface  air  temperature  record 
shown  in  Figure  23.  The  data  have  been  normalized  to  zero  mean  and  unit 
variance.  The  dotted  line  gives  the  expected  value  of  the  range  for  a  normal 
variate. 
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Figure  27.  Histogram  for  the  global  annual  average  surface  air  temperature  record  after 
removal  of  a  trend  determined  by  a  least  squares  fit.  The  residuals  have  been 
normalized  to  unit  variance.  The  dotted  line  is  the  expected  distribution  of  a 
normal  variate. 
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Figure  28.  Variation  of  the  range  of  the  residual  temperature  obtained  by  removing  the  trend 
from  the  data  displayed  in  Figure  23.  The  curve  represents  the  expected  range  for 
a  normally  distributed  variate. 
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a  statistical  sense.  Similar  considerations  hold  for  models.  The  models  may 
generate  series  that  for  time  scales  of  a  1000  years  have  statistics  identical 
to  those  of  a  random  series.  Using  the  models  for  deterministic  prediction 
would  in  this  case  not  be  possible. 
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8  EFFECT  OF  A  SHIFT  OF  THE  MEAN 
ON  THE  FREQUENCY  OF  EXTREMES 


The  consensus  view  is  that  the  greenhouse  effect  will  bring  about  an 
overall  warming  of  the  planet,  together  with  increased  precipitation.  But  it 
is  also  anticipated  that  there  will  be  strong  regional  variations  with  some 
regions  cooling  and  others  receiving  less  precipitation.  The  primary  effect  of 
these  changes  will  be  to  shift  the  mean  although  there  may  also  be  changes 
in  the  variance  other  than  those  resulting  from  trends.  The  shift  in  mean  can 
cause  a  large  change  in  the  probability  of  extreme  events,  as  is  evident  from 
Figure  2.  If  over  100  years  a  particular  climate  parameter  shifts  by  one-half 
a  standard  deviation  toward  lower  or  higher  values,  the  initial  probability 
of  0.5  becomes  0.05  or  0.95,  respectively.  The  large-scale  wind,  temperature 
and  moisture  patterns  in  the  atmosphere  and  the  land  and  ocean  changes 
associated  with  them  undergo  variations  at  all  time  scales  from  hours  to  bil¬ 
lions  of  years.  Some  long-term  changes  have  their  origin  in  events  that  are 
external  to  the  atmosphere-ocean  system,  such  as  variation  in  the  earth’s 
orbital  parameters,  volcanic  explosions  and  changes  in  atmospheric  compo¬ 
sition.  These  long-term  changes  can  accentuate  or  ameliorate  the  impact  of 
shorter-term  fluctuations  which  have  their  origin  in  chaotic  noise.  Alterna¬ 
tively,  long- period  chaotic  noise  might  alter  the  local  mean  and  carrys  the 
short  term  fluctuations  so  that  extreme  values  are  reached.  In  this  section  we 
briefly  consider  the  effect  of  shifts  of  the  mean  on  the  frequency  of  extremes. 


The  events  of  the  Little  Ice  Age  illustrate  the  impact  of  shifts  in  the 


mean  on  extreme  events.  Parry  and  Carter  (1985)  have  examined  in  detail 
the  effect  of  short-term  cool  spells  during  the  Little  Ice  Age  on  marginal 
agriculture  in  southwestern  Scotland.  In  particular,  they  showed  that  small 
changes  in  the  mean  temperature  associated  with  the  Little  Ice  Age  increased 
the  frequency  of  damaging  weather.  Further,  the  probability  of  two  successive 
bad  years  is  even  more  sensitive  to  changes  in  the  mean.  The  importance 
of  successive  extremes  lies  in  their  cumulative  impact:  an  agricultural  region 
may  be  able  to  understand  a  single  shock,  but  if  buffer  stores  are  depleted 
by  one  bad  season,  a  second  one  in  succession  may  be  far  more  devastating. 

The  longest  available  temperature  record  available  to  us  is  the  one  com¬ 
piled  by  Manley  (1974)  for  Central  England  from  several  discontinuous  series. 
The  record  shows  a  number  of  isolated  cool  years  (1740,  1782,  1860,  1979, 
1922)  and  a  clustering  of  two,  three,  or  even  more  extremes  in  successive 
years  (1673-75, 1688-98, 1838-40,  1887-88,  1891-92)  as  indicated  in  Figure  29. 
During  the  clustering  of  several  cool  summers  in  a  row,  when  the  tem¬ 
perature  was  marginal  for  the  growth  of  oats,  the  farms  at  higher  elevations 
in  Scotland  suffered  greatly.  One  cool  period,  which  persisted  from  1688  to 
1698,  led  to  the  abandonment  of  most  of  the  farms  above  300  m  (Parry  and 
Carter,  1985).  Indeed,  the  “Seven  Ice  Years”  of  the  1690s  caused  catastrophe 
among  rural  populations  in  all  of  Northern  Europe. 

The  effect  of  a  change  in  the  mean  on  the  probability  of  an  extreme  event 
can  be  understood  from  a  simple  statistical  model.  The  probability  that 
some  climate  parameter  (e.g.,  summer  temperature,  spring  precipitation) 
falls  above  or  below  some  threshold  value  for  damage  is  P\.  In  a  stable  climate 
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Figure  29.  Variations  in  yearly  average  temperature  in  Central  England.  After  Manley 
(1973). 


65 


regime  the  return  period  for  the  damaging  event  is  Pf 1  years.  Suppose  that 
the  climate  shifts  in  such  a  way  that  the  mean  of  the  climate  parameter 
shifts  by  X  standard  deviations  toward  the  threshold  value  and  that  the 
climate  parameter  is  normally  distributed.  The  ratio  of  the  probability  in 
the  shifted  regime  to  that  in  the  original  regime,  Ps/Pi ,  is  a  highly  non¬ 
linear  function  of  the  shift  of  the  mean,  as  is  shown  in  Figures  30  and  31. 
If  the  initial  probability  is  Pj  =  0.05  (a  return  period  of  20  years)  and  the 
mean  is  shifted  by  one  standard  deviation  toward  the  threshold,  then  the 
subsequent  probability  is  Ps  —  0.25,  or  a  return  period  of  four  years.  The 
probability  of  two  successive  years  of  an  event  with  a  return  period  of  20 
years  is  only  1  in  400  in  the  initial  state  but  moves  to  1  in  16  in  the  shifted 
state.  For  example,  the  standard  deviation  for  the  temperature  record  of 
central  England  is  0.58°C,  with  a  mean  for  the  1659-1977  period  of  9.1°C. 
A  decrease  in  the  mean  to  8.5°C  would  alter  the  probability  of  a  once-in- 
100  years  cool  year  to  a  once- in- 12- years  event.  The  long- period  fluctuation 
in  climate  that  comprised  the  Little  Ice  Age  increased  the  probability  of  a 
short-term  cool  summer,  making  a  sequence  of  successive  cool  summers  much 
more  likely  than  in  the  preceding  or  subsequent  years. 

Changes  in  the  frequency  of  extreme  high  or  low  temperature,  or  of 
high  or  low  precipitation,  can  be  expected  as  climate  shifts  in  response  to 
greenhouse  warming.  The  observed  change  in  global  average  temperature  of 
about  0.5°  over  the  last  century  represents  a  shift  of  2.3  standard  deviations, 
making  the  probability  of  extremely  warm  global  average  years  significantly 
higher. 
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Figure  30.  Ratio  of  probability  of  an  extreme  event  to  the  probability  of  the  event  in  the 
original,  stable  climate.  Shift  in  mean  measure  is  in  terms  of  standard 
deviation.  A  normal  population  is  assumed. 
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Shift  in  Mean 

Figure  31.  Shift  in  probability  from  a  change  in  mean  for  initially  rare  events.  Shift  in 
mean  is  measured  in  terms  of  standard  deviations,  and  a  normal  population 
is  assumed. 


9  SUMMARY 


The  statistical  theory  of  extreme  events  can  be  applied  in  numerous 
ways  to  problems  in  climate  change.  Calculation  of  the  range  as  a  function 
of  time  can  be  used  to  detect  the  presence  of  a  trend.  The  value  of  using 
extreme  events  in  this  application  is  that  since  they  are  rare,  extreme  events 
can  be  treated  as  independent  without  any  assumption  with  respect  to  the 
underlying  distribution.  Analysis  of  the  global  annual  average  temperature 
record  using  extreme  events  and  analysis  based  on  an  assumptions  of  an 
underlying  normal  distribution  both  confirm  the  existence  of  a  trend  with 
high  odds  (10s  —  106)  in  favor  of  the  trend  hypothesis. 

Analysis  of  observed  long-term  global  annual  average  surface  air  tem¬ 
perature  records  and  of  model  results  show  that  both  kinds  of  records  are 
indistinguishable  from  a  series  of  normally  distributed  variates  once  a  trend  is 
removed.  For  both  kinds  of  records,  the  statistics  approach  those  of  normally 
distributed  variates  about  the  mean  and  in  the  tails  of  the  distribution.  For 
model  results,  this  conclusion  holds  even  for  series  that  are  several  thousand 
years  long.  These  observations  cast  doubt  on  the  use  of  General  Circulation 
Models  to  predict  future  climate  deterministically  . 

The  behavior  of  extreme  values  in  nonlinear  systems  can  be  understood 
in  terms  of  the  physical  dimensions  of  the  underlying  attractor.  If  the  govern¬ 
ing  equations  contain  conservation  laws,  the  outer  boundary  of  the  attractor 
lies  within  a  manifold.  The  dimensions  of  the  manifold  define  the  limits  to 
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the  extreme  values. 
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